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We present a canonical mapping transforming physical boson operators into quadratic products of 
cluster composite bosons that preserves matrix elements of operators when a physical constraint is 
enforced. We map the 2D lattice Bose-Hubbard Hamiltonian into 2x2 composite bosons and solve 
it at mean field. The resulting Mott insulator-superfiuid phase diagram reproduces well Quantum 
Monte Carlo results. The Higgs boson behavior along the particle-hole symmetry line is unraveled 
and in remarkable agreement with experiment. Results for the properties of the ground and excited 
states are competitive with other state-of-the-art approaches, but at a fraction of their computational 
cost. The composite boson mapping here introduced can be readily applied to frustrated many-body 
systems where most methodologies face significant hurdles. 
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Introduction.- In the past few years, there has been 
great experimental progress on the control and manipu- 
lation of cold atomic gases loaded in optical lattices, lead- 
ing to quantum simulators of the Bose-Hubbard model 
and its Mott insulator to superfluid transition [T]. A 
notable recent experiment has revealed the Higgs boson 
behavior across this transition in a 2D optical lattice [2] ■ 
There is currently great interest in cold atomic physics 
for engineering synthetic gauge fields that induce topo- 
logical phases and phase transitions. This can be accom- 
plished using a combination of laser-induced tunneling 
with supcrlattice techniques [3] , or by time-periodic shak- 
ing of the lattice [4|. From the theoretical perspective, 
traditional mean field approaches can describe the phase 
diagram of bosonic atoms in lattices of various geome- 
tries, but only qualitatively jSJIS]. Quantum Monte Carlo 
(QMC) yields a highly accurate description of ground 
state properties at zero and finite temperatures when- 
ever the system has no frustration |S] • Static and dy- 
namic properties have also been studied with the Vari- 
ational Cluster Approximation (VCA) [51 [TO]- Exten- 
sions of static mean field approaches involving the use 
of clusters have been considered [TlTll3j . In this work, 
we introduce a theory that maps cluster subspaces of 
the original Fock space onto composite bosons containing 
the exact internal dynamics of the cluster, and whose in- 
teractions account for residual correlations between the 
clusters. Because the mapping is canonical, it is then 
possible to apply standard many-body techniques to this 
Composite Boson (CB) Hamiltonian. In this sense, the 
present model builds upon Hierarchical Mean Field The- 
ory (HMFT) for quantum magnetism [T3]. These ideas 
are here extended to interacting bosons systems loaded in 
optical lattices. We refer to the resulting model as Com- 
posite Boson Mean Field Theory (CBMFT). We demon- 
strate that the inclusion of higher order fluctuation terms 
in the composite mean field yields very accurate results. 
The CB approach to the Bose-Hubbard model unravels 
the Higgs boson behavior along the particle-hole sym- 



metry line and yields remarkable agreement with exper- 
imental data [2]. 

Composite Boson Mapping.- Let us start our deriva- 
tion by assuming a quadratic mapping of the elementary 
boson creation (annihilation) operators a\ (a^) at site i 
in terms of a generic set of CB W a (b a ). 
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If the states |a) of the CB space are known, the matrix 
U is formally defined as U* a/3 — (a\ a\ \0) and Uip a — 
(a\ ai and has unitary character in the Wigner defi- 
nition. Let us now explore the conditions which should 
be fulfilled by transformation ([!]) in order to preserve the 
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canonical bosonic commutation relations a, , a J 

|_ 3 

Inserting the transformation in the commutator we ob- 
tain 



h °5] = E IPiofViaF - Uifl'aU*^) b\bp, = (2) 
a/3/3' 

J2 ({PI ^ \a) (a\ a] \?) - (/?| a] \a) (a\ a t |/3')) b^p,. 



The satisfaction of the canonical commutation relations 
rely on i) resolution of the identity, J2 a \ a ) ( a \ ~ I> anc ^ 
ii) fulfillment of the physical constraint, b\jj a = I. 
Both conditions are equivalent and can be satisfied by 
choosing a complete set of states of the composite Hilbert 
space. However, the mapping has no significant power at 
this stage, as we would still have to deal with the whole 
Fock space of the system. Therefore, we divide it into 
clusters and derive an exact Hamiltonian by means of 
the corresponding cluster composite bosons. Following 
this idea, any operator that is an algebraic function of 
the original operators a, within a single cluster will be 
mapped to a one-body operator of the CB space. In the 
same way, any product of operators belonging to N dif- 
ferent clusters will be mapped to an iV-body operator. 
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Therefore, a general two-body boson Hamiltonian would 
lead to a four-body CB interaction. For the sake of sim- 
plicity we will here restrict ourselves to a density-density 
interaction that leads to a two-body CB Hamiltonian, 
since each density operator is contained in a single clus- 
ter. Let us assume a general Hamiltonian of this kind 



H = 



\tija\cij 



+ VijUiUj 



a] a*. 
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We start our derivation assuming a square lattice par- 
titioned into a set of M clusters, each one at position 
R of a CB superlattice and containing L x L sites. We 
assume a perfect tiling between the superlattice of CB 
bosons and the original lattice. Next, we formally map 
the Hamiltonian using the prescription described above 
and rewrite it in terms of CB labeled by the occupation 
configuration of each cluster, n = {ni, . . . , til, ■ ■ ■ ^l 2 }, 



Hcb = ^2 {Rn\H\Rm)b^ Rn b Rl . 
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* &k & fl'n'W>fl'm'- (4) 



For reasons that will become clear below, we will per- 
form a generic unitary transformation among the CBs, 
b Ra = S n kfln^Rn- I n ^his particular case, we are inter- 
ested in a uniform and real ground-state wavefunction, 
therefore, we will drop the superlattice site index R from 
now on. In this new basis, the hamiltonian can be written 
as 



rti3 



R a/3 



RR' aa'PP' 

where the matrix elements 



E ^Urf^^. (5) 
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E (nn'|ff|mm') COT/ (7) 



1 1 1 1 i r 1 1 1 1 



encode all the information of the original Hamiltonian. 
Notice that, in general, (nn'| H |mm') depends on the 
precise layout of the two considered neighboring clusters. 
The CB Hamiltonian ^ is an exact image of the original 
boson Hamiltonian provided that the physical constraint 
in each cluster site (J2 a "RcfiRa — I) is satisfied. On the 
other hand, treating this Hamiltonian by means of stan- 
dard many-body techniques, we immediately incorporate 
quantum correlations inside the cluster in an exact way. 

Composite Boson Mean Field Theory - We here treat 
the CB Hamiltonian in the Bogoliubov-de Gennes ap- 
proximation. In order to proceed further, we have to 



specify the matrix elements of the initial lattice Hamilto- 
nian. As a first test of CBMFT, we benchmark the Bose- 
Hubbard Hamiltonian in a 2D square lattice. Namely, 
tij = —t5i j +e where e is the unit vector in the lattice 
directions x, y, and = USij is the on-site Hub- 
bard repulsion. In what follows, we omit U and mea- 
sure all quantities in units of U. Assuming a uniform 
2D lattice with translational symmetry, we first per- 
form a Fourier transform of the CB boson operators 



(l/VA/)£ k e~ 



iik-R tt 



Ka' leading to 



H cb = E^XXcA/s (8) 

a/3 k 

+T7 E W W' E Tq b L* b b+<i«' b ki+ii3 b W> 



aa'00' 



k x k 2 q 



where we have introduced 7 q = cos(Lq x ) + cos{Lq y ) after 
a symmetrization of the two-body matrix elements W^p 
in order to preserve the lattice C4 symmetry. Details on 
the calculation of these matrix elements can be found in 
the Supplemental Material. Next, we assume a conden- 
sation of the CB in the k = 0,a — state by introduc- 
ing a shift transformation &k=o,a=o = fyk=o q=o = <J V / M. 
This transformation manifestly violates the on-site physi- 
cal constraint. Thus, we relax it and impose a global con- 
straint on the CB density, X^Sa b R a b Ra = M. Trans- 
forming to momentum space and shifting, the global 
physical constraint becomes 



OVr, = 1 



ka°ka 



(9) 



a^O k 



where we have neglected the fluctuations of the con- 
densed boson. Eq. (JoJ) defines a 2 as the CB conden- 
sate fraction. Inserting the constraint ([9]) by means of a 
Lagrange multiplier A in the CB Hamiltonian ([8]) and ap- 
plying a mean field decoupling, we arrive to a quadratic 
Hamiltonian of the form 



H MF = + J2 A ^b{ a b kl3 

ka/3 

+ E ( B ^ b l a b -kp + #ka/?&-k/?&ka) (10) 



kail 



where the specific form of H^ ' and the matrices A ka a 
and i?ka/3 can be found in the Supplemental Material. 
The quadratic mean field Hamiltonian ( 10 1 can be 



diagonalized by means of a Bogoliubov transformation 
"L = S,s^k,Sa&k/3 _ llpYhab-kp, leading to the Bo- 



goliubov eigenvalue equation [15] 
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where the positive eigenvalues uj k determine the excita- 
tion spectrum. The Bogoliubov equation depends on the 
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FIG. 1: (color online) Phase diagram of the first lobe of the 
Mott-insulator to superfiuid transition. The dotted, solid, 
and dashed curves show the composite 2x2 boson mean field 
results in increasing order of approximation. The particle- 
hole symmetry line traverses the Mott and superfiuid regions. 
Black circles are QM C results from Ref . 17\ . Grey circles dis- 
play three points in parameter space for which the dispersions 
are analyzed in Fig. 3. 



generic transformation U previously defined. We are now 
ready to derive a Hartree-Bose (HB) like equation for this 
transformation upon minimization of the free energy with 
respect to the condensed CB structure U^. The resulting 
equation can be cast in matrix form 



(12) 



The derivation of the matrix elements of h for L x L 
clusters, given in the Supplemental Material, is straight- 
forward though lengthy. The HB Hamiltonian h depends 
on the Hartree transformation U , on the Bogoliubov am- 
plitudes Xk, lit, and on the fraction of the conden- 
sate a 2 . Strictly speaking, the self-consistent Hartree 
diagonalization provides a single eigenvector defining the 
structure of the condensed CB and the corresponding 
lowest eigenvalue, which is the Lagrange multiplier A. 
However, after attaining self-consistency the matrix di- 
agonalization procedure supplies a complete set of eigen- 
states that are orthogonal to the condensed CB. It is in 
this basis orthogonal to the condensate where the Bo- 



goliubov Hamiltonian ( 10 1 is expressed. We seek a self- 
consistent solution of the coupled set of equations de- 
termined by the diagonalization of the HB Hamiltonian 
(12 1 that fixes the Hartree transformation U and the 



Langrange multiplier A, and the Bogoliubov equations 



(111 that provides the Bogoliubov amplitudes X^ and 
Yk, together with the expectation value of the physical 
condition ^ for the determination of the CB condensed 
fraction a. 



2D Bose-Hubbard Model Results- We start with 
benchmark calculations describing the insulating Mott 
phase in the first lobe corresponding to a density per site 
of p = 1 and cluster size of 2 x 2. Within this phase, the 
structure of the Hartree transformation U and the CB 
condensed fraction a 2 are //-independent. The structure 
of the condensed CB dictated by U° is a linear combina- 
tion of cluster states with N = 4 physical bosons. The 
CB fluctuations are pairs of particle (N — 5) and hole 
(N = 3) physical bosonic states. In addition, particle- 
and hole-like excitation eigenvalues have a linear depen- 
dence on the chemical potential for fixed t. Both excita- 
tions cross each other at the particle-hole symmetry line, 
where the gap is doubly degenerate. The edges of the first 
Mott lobe are determined by the vanishing of the gap, in- 
dicating the appearance of a Goldstone mode at k = 
related to the U(l) symmetry breaking in the superfiuid. 
Fig. 1 shows the phase diagram of the Bose-Hubbard 
model in three different CB mean field approximations. 
The th order approximation neglects fluctuations and 
solves the Hartree equations exclusively. The edges of the 
Mott lobe are determined in this case by a deviation from 
the density p = 1. This th order approximation is equiv- 
alent to the cluster mean field calculations of Refs. [TlT 
H3| . The 2 nd and i th order approximations incorporate 
fluctuations by means of a self-consistent solution of the 



Bogoliubov (11) plus Hartree ( 12 ) equations linked by the 
physical condition The 2" d order approximation ne- 
glects two-body interactions between fluctuating bosons, 
while the 4 th order solves the three coupled equations 
in full. As the approximation order increases, CBMFT 
shows clear convergence towards QMC. Also shown in 
Fig. 1 is the extension of the particle-hole line to the 
superfiuid phase characterized by density p = 1. 

The full self-consistent A th order approximation does 
not describe the gapless feature of the superfiuid phase 
correctly. Although ways to correct this deficiency have 
been suggested [TB] , in the rest of this paper we will focus 
on the 2 nd order approximation that strictly preserves the 
gapless spectrum when U(l) symmetry is broken. 

Fig. 2 shows the total density p = 

condensate density p c — |(0|a||</>)| 2 for hopping values of 
t = 0.02 and 0.04. VCA results for t = 0.02 [HUTU] com- 
pare well with our results in Fig. 2. The plateau charac- 
terizing the Mott phase is reduced for larger t. Outside 
this region, the superfiuid has non-commensurate density 
in the entire phase diagram except for the particle-hole 
(p-h) line. The condensate density of physical bosons, 
representing the coherence of the superfiuid phase, van- 
ishes in the Mott phase. 

In Fig. 1, we have depicted three characteristic points 
at t = 0.04; namely, a is at the p-h line in the Mott phase, 
b is still in the Mott phase but away from the p-h line, 
and c is in the superfiuid phase. Fig. 3 shows particle- 
and hole-like excitations for k v = as a function of k x for 
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FIG. 2: (color online). Total density (p) and condensate den- 
sity (p c ) for t — 0.02 (solid line) and t = 0.04 (dashed line) 
within the second order approximation. 



FIG. 4: (color online). Higgs, Goldstone, particle, and hole 
modes along the p-h symmetry line computed within 2 nd or- 
der CBMFT (solid line). Experimental data points from |2]. 
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FIG. 3: (color online). Dispersion modes for t = 0.04: 
particle- and hole-like excitation modes at (a) the particle- 
hole symmetry line(/i = 0.419), and (b) inside the Mott insu- 
lator (n = 0.30), and amplitude- and phase-like modes at (c) 
in the superfluid (/x = 0.12). 



a and b inside the Mott phase. The degeneracy of the 
particle and hole modes for point a approaching k = 
is clearly seen in this figure. Away from the p-h line 
and still in the Mott phase (point b), this degeneracy is 
broken and the hole is favored against the particle mode. 
Well inside the superfluid phase (point c) , we recognize a 
gapless mode (Goldstone) with the characteristic linear 
dispersion at low momentum, as well as a gapped mode. 
An analysis of the CB structure of each mode, similar 
to the one performed in Ref. [TT] shows that the gapless 
mode is a phase-like mode, while the gapped mode is an 
amplitude-like mode. 



As already shown in Fig. 1, the p-h line crosses the 
tip of the Mott lobe where the quantum phase transi- 
tion is driven at constant density (p = 1). The phase 
transition at the lobe tip can be understood within the 
Higgs mechanism, as has been recently discussed in Refs. 
[21 [17] . Fig. [1] displays how doubly degenerate excita- 
tions along the p-h line inside the Mott insulator vanish 
at the critical point. In the superfluid region, one of them 
remains at zero excitation energy (Goldstone) while the 
other one grows for increasing hopping (Higgs). In both 
cases, their structure mixes particle- and hole-like states 
of the cluster. The CBMFT results not only match the 
experimental data [5] remarkably well but also gives an 
excellent description of the critical point. 

Conclusions.- We have introduced a cluster composite 
boson mapping which separates intra- and inter-cluster 
degrees of freedom. The former are treated exactly while 
the latter can be approximated using standard many- 
body methods applied to the resulting CB Hamiltonian. 
We have here shown that a mean field approximation to 
the CB interaction for the Bose- Hubbard model produces 
an accurate description of the Mott-superfluid phase di- 
agram compared to QMC results. Densities and disper- 
sions are found in quantitative agreement with more so- 
phisticated techniques like VCA. The recently measured 
Higgs mode is also computed and found to be in remark- 
able agreement with experiment. Further improvement 
of the theory beyond the mean field second-order ap- 
proximation employed in this paper is feasible. Most im- 
portantly, CBMFT is readily applicable to other many- 
body problems where frustration, synthetic gauge fields 
or long range interactions pose significant hurdles to ex- 
isting state-of-the-art methodologies. 
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SUPPLEMENTAL MATERIAL 



Matrix Elements of the CB Hamiltonian 



The cluster-cluster interaction matrix elements depend on the particular layout of the two clusters R and R' . For 
the Bose-Hubbard model, it is exclusively related to the hopping term that destroys (creates) a boson at the edge of 
one cluster and creates (destroys) a boson at the edge of a neighbor cluster. The specific form of the interaction for 
L x L clusters is 



<i,j) n,n' 



+ \/rii\ / n\ + 



1 U R{ni-l} U R> {n^+ljU R n U R'n' 



(13) 



where Yluj) rvms over the L links between the edges of both clusters with i g R and j € R' . Using compact notation, 
{rii + 1} corresponds to a configuration n with one more boson at site i, that is, {rij + 1} = (n 1; . . . , + 1, . . . , n L i). 



Next, we perform the C4 symmetrization of the matrix elements (131 in order to preserve the lattice symmetry at 
the mean field level 



W $ = \ E E (v^+iy^ u tn i +i} u ti- j -i} u ^ + Vm^j + 1 ^Vi}^ + i } ^5) , (14) 

n,n' 

where now Y]uj\ runs over the L x L links which connects cluster R with the four nearest neighbors R' . The 
intra-cluster one-body matrix elements has both contributions, the hopping term and the Hubbard interaction 



T t = E E (h fa - J ) - ^) u ^ 



n jeR 



-*E E (V^+iV^ Ul ni+lin ._ 1} u£ + y^v^TTT *7 { Vx, n . +11 £^ 

n (ij)eR 



(15) 
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Bogoliubov matrix 

The components of the mean field Hamiltonian of Eq. (10) in the Letter, H^ ' and the matrices A ka p and Bk a p, 
are given by 

H<® = 2MW ° V + M (T ° - A) a 2 

A ka0 = (T| - XS aP ) + 2a 2 (Wggn + 2W$) + ^ E E (^' r kq + 2W^) P qa ^ 

B* afj = a 2 W^ lk + 1 E E <^ kq iW (16) 

q a'/3' 

where we have defined Tk q = cos(Lk x ) cos(Lq x ) + cos(Lk y ) cos(Lq y ). The density matrix P qa( g = (& qa /fr q /3' ) = P q /3 a 



J qa' u <iP I ~ ± qfi 

is hermitian, and the pairing tensor Kq a p — (6q a /&l qi g/^ = ^- q /3a is symmetric. For the particular case of the Bose- 
Hubbard model, we assume that both of them are real. 

Both matrices A k and B k depend on the density matrix and the pairing tensor, which were defined as expectation 
values in the Bogoliubov quasi-particle vacuum. Upon inversion of the Bogoliubov transformation, it is straightforward 
to show that 

Pka/3 = E ^kaa'^k/3a', K ka p = ^ Y\. aol i X k fj a t . (17) 



Hartree-Bose matrix 

Minimization of the free energy with respect to the condensed CB structure (see Eq. 18 in the letter), leads 
to the Hartree-Bose (HB) eigensystem: 



E 



hm, n U° = E (h£]n + = (18) 



where we have separated the zeroth and second order contributions to the HB Hamiltonian. 

The zeroth order has no contribution from fluctuating bosons, while the second order has quadratic contributions 
from fluctuating bosons 

(JT<°>) = 2MW ° V + M (T ° - A) a 2 , (19) 
(H^) = 2a 2 (wSf-ykKtafi + [W$ 7 u + 2W$] P ka/3 ) ■ (20) 



Making use of Eqs. ( 14 1 and ( 15 ) we first write explicitly the form of the tensors for the interaction and self-energy 



of the condensed CB entering in Eq. ( 19 ) 



V n / \ n' / 



n jeR 



- f E E (^H^^ +1 ,„ ^1) ^/^^v , ^li-l,n J+ l)4 (22) 

n {ij)£R 

The corresponding derivatives with respect to the condensate CB structure function are 
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i Y, (V^Ti u° {mi+1} + u° {nH _ 1} ) [J2 + 1 u {n- j+ i}Ut> J 

<ij) V n' / 

+ ^ E ( E^^ 1 ^{n i+ l}^n J (^TTT f/ { ° ro . +1} + ^ E/° m ._ 1} 

E ( E ^K+i}^ J (y^TTT t/ { ° m . +1} + ^ U° {mj _ 1} ) , 



(23) 



and 



T = 2 E ( ~ X ) ~ ^ ) ^m.n 



Collecting terms, the zeroth order HB matrix is 



-2t Y (y/™iy/ m j + 1 U { mi -i,m j+ i} + + ly/m] Uf m . +l m;j _ 1} ) . 



(24) 



^min = E ( ^3 ~ ^ ~ ^ ) 5m > a 



-t Y (V^iV m 3 + 1 (5 n,{m,-l,m J + l} + V™* + 1\A™J ^,{m, + Lm r l)) 



(»3>6-R 



(ij) L V n' / 



+ 1} + V™>j ^)m r l}) 



(25) 



Following a similar procedure, we obtain the second order contributions to the HB matrix (20). First we write 



explicitly the tensors for the interaction between two condensed bosons and two fluctuating bosons, and then their 
corresponding derivatives 



a/3 



III) 



w a0 



vv 0/3 — 



4 E E i^+^S u? 7H+1} u? n ,_ 1} u°K, + u? ni _ iy u? n , +1} uZK,) , (26) 

(ij) n,n' 

\ E E (^7+1^7 U? ni+1} U? m ._ 1} + ^l^+l t/f rii _ 1} C/f mj+1} + [n o m]) Ul (27) 



BUI 



\ E E (%A^ + T^ + 1 f/f ni+1} C/°,C/°C/f n , +1} + U? ni _ 1} U°,U°U? n ,_ 1} ) , (28) 

\ J2 E + ^{« <+ i} + V^V^ Ul ni _ 1} U? mj _ 1} + [n o m]) (29) 



'+l}^n' J ; 



1 



(ij) \ n' 

Making use of the W tensor derivatives, the second order term of the HB matrix is 



(30) 
(31) 



(ij) 

+ E (V^+Ty/mj + 1 U? ni+1} U* mj+1} + V^iV^J U? ni _ 1} U? mj _ 1} + [n^ m]) 

<ij> 



k 

X] 7k-Pka/3 



(32) 



+ 2AfE 



J L k 



